Extrema statistics in the dynamics of a non-Gaussian random field 
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When the equations that govern the dynamics of a random field are nonlinear, the field can 
develop with time non-Gaussian statistics even if its initial condition is Gaussian. Here, we provide 
a general framework for calculating the effect of the underlying nonlinear dynamics on the relative 
densities of maxima and minima of the field. Using this simple geometrical probe, we can identify 
the size of the non-Gaussian contributions in the random field, or alternatively the magnitude of the 
nonlinear terms in the underlying equations of motion. We demonstrate our approach by applying 
it to an initially Gaussian field that evolves according to the deterministic KPZ equation, which 
models surface growth and shock dynamics. 



Random fields that undergo a time evolution accord- 
ing to a nonlinear dynamical equation often develop non- 
Gaussian statistics that provide clues about the details of 
the underlying microscopic mechanisms. Consider for ex- 
ample a gas-liquid phase transition. In the early stages, 
there are many randomly small volumes in which all the 
molecules are in the same phase, distributed randomly. 
Over time, these volumes will grow and merge, thereby 
gadually replacing the Gaussian disorder with structure 

Even if the initial condition of a random field is Gaus- 
sian, the dynamics will typically generate a non-Gaussian 
component in the field that we wish to quantify and track 
with time. The standard approach to detect and measure 
non-Gaussianities is to employ higher-order correlation 
functions. In this work, we adopt a geometric approach 
to measuring the non- Gaussian component of a scalar 
field h(r,t): we interpret it as a height function describ- 
ing an evolving surface, and study its geometry. Gaussian 
surfaces have certain general geometric and topological 
properties 043- For example, the number of maxima ex- 
actly balances the number of minima. A random surface 
that does not exhibit this property is then guaranteed to 
have non-Gaussian statistics 



perty 



In previous articles 0, H| we studied fields that arc 
local functions of a given Gaussian, i.e. of the form 
h(r) = H(r) + f n l{H (r)) , where H is a Gaussian field 
and Jnl a nonlinear function. In this scheme, the per- 
turbed height h at any point r is a function only of the 
original height H(r) at the same point. In this paper, we 
move to the general case of nonlocal perturbations, which 
e.g. include a dependence on VH, thereby introducing a 
mixing between the field values at different points. 

Such a nonlocal non-Gaussianity can arise in a broad 
range of physical contexts, for example as the result of 
nonlinear diffusion. For concreteness, consider a diffusion 



equation of the general form 



dh(r, t) 
dt 



= DV 2 h(r 1 t) + f NL {h,Vh), 



(1) 



where Jnl is any nonlinear function. If we let h be 
a Gaussian field at t = 0, then non-Gaussianities will 
emerge as a consequence of the last term; if we would 
omit this term, we retrieve the heat equation, which 
would preserve the Gaussianity of h for all t > 0. A 
variety of known diffusion equations has this general 
form. For instance, when /jvl takes the form — h 2 we 
get Fisher's equation, which can be used as a model to 
describe the growth and saturation of a population. An- 
other example is the Cahn-Hilliard equation for the de- 
velopment of order after a phase transition [![. Several 
models of structure formation, in both condensed matter 
@ and cosmology [loj], also belong to this class. 
To illustrate our general result, we apply it to the case 



of a field obeying the deterministic KPZ equation [111 ], 
for which /jyx = -|(V/i) 2 . This equation is often used to 
model the height profile of a growing surface. A field that 
starts out as a Gaussian field will acquire non-Gaussian 
characteristics as time progresses. We use our formula 
to quantify the resulting effect on the relative difference 
in densities of maxima and minima. This allows to back 
up the non-Gaussian component in h, or alternatively, to 
deduce what the nonlinear coefficient A is. We verify the 
analytical predictions by comparing them with results 
from computer simulations. 

The outline of this paper is as follows. In section U 
we determine a general expression for the imbalance be- 
tween maxima and minima for a non- Gaussian field. This 
is applied to the KPZ equation in section [TTJ Finally, sec- 
tion [TTT] summarizes our conclusions. 



I. NON-GAUSSIAN FIELDS 

A homogeneous and isotropic Gaussian field is defined 
in terms of its Fourier components as 
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H(f) = ^2A(k)coB(k-f +</>%). (2) 
% 
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The phases <f>^ are independent random variables, uni- 
formly distributed between and 2tt. The amplitude 
spectrum A(k) depends only on the magnitude of the 
wave vector k and encodes the special features of the 
Gaussian field under consideration. An alternative ap- 
proach is to express the amplitude spectrum in terms of 
its moments, according to 

K n = h A (k) 2 k n . (3) 

k 

For convenience, we will consider H to be normalized, 
such that K a = (H 2 ) = 1, see ref. Q for more details. 

In what follows, we concentrate on homogeneous and 
isotropic fields h(f), which we assume to be in the form 
of a Gaussian H[f) with the addition of a perturbation. 
Unlike refs. @, H|, we will not restrict ourselves to a per- 
turbation of the local kind, i.e. where the perturbation 
at any point r is a function of H(f) only. We will now 
also accommodate perturbations which depend on V H 
for instance, or evolve over time. Such perturbations in- 
troduce a mixing between the values of the field at dif- 
ferent points, which we will designate as nonlocal pertur- 
bations. 

We will investigate the effect of a perturbation on the 
densities of maxima and minima. A maximum (min- 
imum) of h is defined by the condition h x {fa) = 
h y (fo) — 0, along with the inequalities h xx (fo)h yy (fa) — 
h xy (fh) 2 > (if this were negative, fo would be a sad- 
dle point) and h xx {fo), h yy (fo) negative (positive); note 
that the first condition implies that h xx (fo) and h yy (fo) 
have the same sign. The x and y subscripts indicate 
derivatives with respect to the coordinates of the two- 
dimensional plane. 

The general procedure that we use is very similar to 
the one in Q and is as follows: we consider a fixed point 
fa - due to the homogeneity of h, the analysis will not 
depend on this choice. We determine the joint probabil- 
ity distribution of h x , h y , h xx , h yy and h xy , since these 
stochastic variables are the ingredients from which max- 
ima and minima arc defined, as outlined above. This 
distribution can be determined via the generating func- 
tion, which in turn can be constructed by determining 
the relevant cumulants involving the five stochastic vari- 
ables. Once the probability distribution is obtained, we 
set h x = h y = and integrate the second derivatives 
over the region defining a minimum (maximum) in order 
to get the density of minima (maxima). 

As we did in [8(, we transform to another coordinate 
system, based on the complex coordinates z = x + iy 
and z* , which will allow us to make full use of the homo- 
geneity and isotropy of h later on. In this new basis, we 
have 

d _ i d d _ l d + ij ^ 

dz dx dy dz* dx dy 

In this coordinate system, the definition of a maximum 
(minimum) becomes h z (fo) = 0, \h zz *(fh)\ > \h zz (fh)\ 
and h zz -{fo) is negative (positive). [l4[ 



Some care is required however, since we are now deal- 
ing with complex variables h z and h zz {h zz * is real). We 
will treat the variables z and z* as if they were inde- 
pendent. Therefore, next to h Zl we will consider h z * as 
well, as a separate random variable, although it is ac- 
tually the complex conjugate of h z . Similarly, we also 
include h z * z * = h* z . Therefore, we are still dealing with 
five variables: h z . h zz , their conjugates, and h zz *. 

As stated before, we will arrive at the joint probability 
distribution of these variables by building the generating 
function, which is the Fourier transform of the probability 
distribution. For a set of n correlated variables £j this is 

x(Ai,...,A„) 

= /di,...,^, <:„),■ ^ •. (5) 

By expanding the exponential into a Taylor series we find 
that the coefficients - which are called the moments of 
the distribution (not to be confused with the moments 
from cq. ([5])) - are correlations: 

x(Ai,...,A„) 

= l + i£&>A i + ^£& 1 &)A il A ia 

+ 3? zJ Zh & 3 > -V' X J2 -V. . + • • • (6) 

jl,j2,j3 

If we do the same for the logarithm of x, we obtain the 
cumulants: 

i 2 

log X = i °i (6 ) X J + 2j ° 2 (&i ' & ) A '; A .- ' 

i 3 

+ 7! C3(£ji , £72 7 €js)^h x h X 3s + • ■ • ( 7 ) 

ji ,32 ,h 

From eqs. ((SJ) and ([7]) it can be derived that the cumu- 
lants can be factorized into moments, for example 

£3(6,6,6) - (666) - (6X66) - (6)(66) , , 
-(6)(66) + 2(6)(6)(6)- U 

If all the cumulants are known, one can reconstruct the 
generating function and from that obtain the probability 
distribution via an inverse Fourier transformation. 

The defining characteristic of Gaussian variables is 
that all cumulants are zero, apart from the second order 
ones (C2). If h were a Gaussian field, then this would 
apply to p(h Zl hzzj h zz *), since the derivatives of a Gaus- 
sian field are themselves also Gaussian fields. Since h is 
non-Gaussian, this is not the case. The first-order cumu- 
lants are still zero; for instance, we have C(h z ) = (h z ) = 
d z (h) = since (h) is constant due to the homogeneity of 
h. The third-order cumulants are however nonzero. We 
will include these and see how they influence the probabil- 
ity distribution and the densities of maxima and minima. 
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In principle, there are infinitely many nonzero cumu- 
lants. However, a field that is generated by a nonlinear 
differential equation, like cq. (JTJ) , typically has small cu- 
mulants of high order. In particular, if /jvl is a quadratic 
function and the initial conditions are Gaussian, then the 
n-th order cumulants scale like f^J 2 (for n > 2) - see ap- 
pendix [3] Therefore we will only need to determine cu- 
mulants up to third order to get the correction to leading 
order. 

The usefulness of the complex variables z and z* be- 
comes apparent when we look for all nonzero cumulants 
of second and third order involving the five variables we 
have. Since h is isotropic, a moment like (h z *h zz ) should 
not change when we rotate the field by an arbitrary angle 
a. Such a rotation would give z — > e la z and z* — > e~ la z*. 
Incorporating these in the derivatives causes the afore- 
mentioned moment to pick up a factor e la . Since we 
argued that the moment should not be affected by the 
rotation, it must be zero. In general, any moment in- 
volving a different number of z and z* derivatives is zero 
by this argument. Since cumulants can be decomposed 
into moments, as depicted in eq. ((8]), the same applies to 
cumulants. 

Furthermore, translational symmetry implies some re- 
lations between the cumulants. From translational in- 
variance it follows that any correlation should be con- 
stant with respect to r. For instance, using the product 
rule, we have 



= d z *(hlh z *) 



2(h z h z *h z 



(9) 



which gives us the relation present in eq. (|10c[) . 

Therefore, there are only a few independent cumulants 
that are (potentially) nonzero: 



v=(\h z \ 2 ), 



p = (\hl\h zz *) = -Uhlh z * z .) 



a=(\h zz \ 2 ) = (h 2 zz ,), 

ft; 

6 = {\h„\ 2 h„.). 



(10a) 
(10b) 



-\{h 2 z ,h zz ), (10c) 



(lOd) 
(lOe) 



In these definitions, the cumulants have been expanded 
into moments in accordance with eq. (|8]); since the first- 
order correlations are zero, as noted before, only the 
third-order correlations remain. We also introduced the 
shorthand notation \h z \ 2 = h z h* z = h z h z * and similarly 
for \h zz \ 2 . Note also that the third-order cumulants, j3, 
7 and S are close to zero when h is close to being Gaus- 
sian, which we assume. On the other hand, a and a are 
nonzero in general. 

We can now construct the logarithm of the generating 
function as prescribed by eq. (0, 



logX = — cr|A 2 | -a|A zz | 

-m\ z \ 2 \ zz * +0(x 2 z x z . 

-i 7 \l zt -i8\\ zz \ 2 \ zz *. 



a 2 ,a zz ) (ii) 



Note that some cumulants appear multiple times in 
cq. ([7]) since the A's can be permuted (if they are not all 
the same); this explains why for instance the term A zz * 
has a prefactor i/6 whereas the prefactor of |A ZZ | 2 A ZZ * = 
A ZZ A Z * Z *A ZZ * is i (due to the 6 distinct permutations of 
the A's). 

We see that x features an exponential of a third-degree 
polynomial, making the inverse Fourier transform - to be 
performed in order to get the probability distribution - 
nontrivial. Remember however that the cubic terms are 
small owing to the near-Gaussianity of h, allowing us to 
make the expansion 



X = [l - */?|A z | 2 A zz . + i/3(\ 2 z \ z * z * + A 2 , A zz ) 

-f 7 AL* -^|A ZZ | 2 A ZZ ." 
x exp(-cr|A z | 2 -a|A zz | 2 - \a\ 2 zz ,). 
The inverse Fourier transform of this gives (l5| 

p(h z ,h zz ,h zz .) 



(12) 



1 + -^h zz , {\h z \ 2 -a)- ^—{h 2 z h z . z , + h 2 z ,h zz ) 
+7r^(h 3 zz , - 3ah zz *) + -^h zz -*(\h zz \ 2 - a) 



1 



TT 2 \/2Traa 3 / 2 



|7iz| 2 /rr-\h zz \ 2 /a-h 2 z!t , /2a 



(13) 



Now that the joint probability distribution of the rel- 
evant derivatives is obtained, we can set h z = h z * = 
this condition defines a critical point. The joint probabil- 
ity distribution measures how likely it is that h z and h z * 
are close to zero for a certain point r. What is needed 
however is for h z and h z * to be exactly zero for a point 
close to r, since we are looking for a density with respect 
to the (a;,y)-plane. For this, we need to go from a prob- 
ability density with respect to h z and h z * to one with 
respect to z and z* (representing x and y). This is ac- 
complished by multiplying p with the following Jacobian: 



J = 



d{h z ,h z ») 



d{z,z*) 



\h. 



(14) 



Now we are ready to set h z = h z * = and integrate p J 
over h zz and h zz * . The range is determined by the type 
of critical point of interest; focus on the minima first. 
For these we must have \h zz \ < \h zz * | and h zz , > 0. The 
integration over h zz is done by integrating over its real 
and imaginary part. Since the integrand depends only on 
the modulus of h zz , we move to polar coordinates. Let 
us define r = \h zz \ and s = h zz *. The integration range 
is then < r < s, and with eq. t|l 3|) we get 

1 



7r 2 \T2/Kaa i / 2 

'Ids 2Trrdr{s 2 -r 2 )e' r2 /°" s2 / 2a 
'o Jo 

S_ 

7'- 



(15) 



x 



aa 6a 3 



1 - — s + -rMs 3 - 3as) + -^rs(r 2 - a) 



(V 
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This integration is pretty straightforward: although the 
range of r is finite, the integrand is a Gaussian multiplied 
by a polynomial that has only odd degrees of r, hence 
it does not give rise to error functions. The resulting 
integral over s is also standard. The final result reads 



a 



2V3ti 



(16) 



For a Gaussian field, we would have /? = 7 = 5 = 0, 
a = and a = jj-Ki- This would give us n m i n = 

K4/(8y/3TTK 2 ), exactly as given in 

To get the density of maxima, the same integrand as 
in eq. (fT5j) needs to be integrated over the range s < 
and < r < — s. However, note that if we make the 
transformation s — > — s, the range of integration is the 
same as in eq. (fPo)) . Furthermore, note that the transfor- 
mation s — > — s in the integrand is equivalent to (3 — > —/3, 
7 — > — 7 and S — > — 5. With this insight, we easily find 
that the expression for n max is the same as the above, 
except with a plus in place of the first minus. 

With this result, the imbalance between maxima and 
minima is found to be 



An = 




(17) 



This is the main result of this paper. As an illustration, 
we shall now use this result to understand the evolution 
of maxima and minima in the context of a differential 
equation describing surface growth. 



II. KPZ EQUATION 

The deterministic Kardar-Parisi-Zhang (KPZ) equa- 
tion [ll[ is given by 



dh 



= v\J l h 



A 



(18) 



This equation is often used to describe the height profile 
of a growing surface: the first term on the right-hand side 
describes the diffusion of particles along the surface, while 
the second term accounts for the assumption that the 
growth is perpendicular to the slope of the surface, while 
h describes the height along the universal up direction 
[13 . This leads to (see fig. [T} 

^ = A V /TT7W=A + ^(V/ 1 ) 2 + ... (19) 

The leading term A is ignored since it is just a constant 
that docs not affect the profile of the surface. 

Another interpretation of eq. (fT8|) is obtained by taking 
the gradient on both sides, which yields 



dv 



vV 2 v + XvVv, 



(20) 




Figure 1: A geometrical interpretation of the KPZ equation 
applied to a growing surface. The surface is assumed to grow 
perpendicularly at a constant rate A. Measured vertically, the 
growth rate is dh/dt « A(l + ±(V/i) 2 ). 



where v = V/i is a velocity field. This is a vector Burger's 
equation which arises in fluid mechanics. The maxima 
and minima of h correspond to sources and sinks of v. 

We will take h(r,t) to be a Gaussian field at t = 0, 
and use our result eq. (|17|) to determine how the non- 
Gaussianitics, which arise and evolve due to the KPZ 
equation, influence the densities of maxima and minima. 

First note that if we would set A = in eq. (fTg|) . we 
retrieve the heat equation, which preserves the Gaussian- 
ity of a field: if we enter h(f, t = 0) = H(r), where H{f) 
is a Gaussian field as given by eq. ([2]), we find that the 
solution is 



h(f,t) 



k 



A(k)e~ k vt cos{k-r + ct) k ). 



(21) 



We find that the amplitudes pick up a factor exp(— k 2 vt) : 
but the phases remain independent. Therefore, even 
though its amplitude spectrum changes, h(f,t) remains 
Gaussian at any time t and the density of maxima and 
minima remains the same, since this is a general property 
of Gaussian fields. 

If we have A ^ 0, h(f,t) no longer remains Gaussian. 
In fact, as we will see, the density of maxima and minima 
is no longer the same. We shall assume A to be small in 
comparison with v, and find out how these densities differ 
as a function of time, using cq. (1171) . For this, we need 
to determine the two- and three-point correlations a, a, 
P, 7 and (5. 

First, we substitute u = exp((A/2i/)/i). Note that, 
since this is a monotonically increasing function of h, the 
maxima and minima of u are exactly the same points as 
those of h. In terms of u, the KPZ equation becomes: 



du 2 



(22) 
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which is simply the heat equation. However, u(r,t = 
0) = cxp((A/2z/)/io) is now n °t a Gaussian field. If we 
assume that A <C v, we have: 



A 



A 



u 



l + ^o + ^-o^ + 0((A/z/) d ). 



2v 



8;/ 



(23) 



Since the leading term, equal to one, has no influence on 
either the maxima and minima or the heat equation, we 
can ignore it. The same applies to the prefactor A of 
the second term. Hence we make a final transformation 



2/' 



(u-1), 



v = h + —ti 2 + O((\/v) 2 ). 

4:1/ 



(24) 



(25) 



Note that v still obeys the heat equation and also shares 
the same maxima and minima with h and u. Moreover, 
we now have v(r, t = 0) in the desired form of a Gaussian 
ho plus a perturbation. Since v obeys the heat equation, 
we can use the corresponding Green's function to write 
down the general solution 

v(r, t) = J d 2 r G(r, r, t)vo(f) 

= I ^^e-^^ + ^ir?), (26) 

where vo(r) = v(r,t = 0). 

We can now calculate the five correlations needed to 
determine An. We will demonstrate the procedure using 
o = {vzifj^Vz'^jt)) as an example. 

a = (v z (r,t)v z ,(r,t)) 

d 2 fid 2 f2d xl G(ri,fi,f) 

S#!G(ra,f2,*)(«b(fi)«o(f2)> • (27) 

r\=r 2 =r 

The brackets represent averaging over all </>g that define 
vq, while the spatial derivatives act only on the respective 
Green's function. The latter gives 



d Zl G(ri,r!,t) = d Zl 



1 (n-fif 
e ivt 



4-Kvt 



1 



((xi - xi) - i(y - yi))e " 



(28) 



n{4vt) 2 

The moment present in eq. (|27l) is 

{vo{ri)v (f 2 )) 



= (/io(ri)/i (f 2 )) 

+ ±((h (r 1 )h (r 2 f) + (Mri) 2 fco(f 2 )» 



+ (^Y(h (h) 2 h (f 2 ) 2 ) 



\4v 



(29) 



Note that the second term (the one linear in \/4v) is 
a three-point correlation, and therefore zero due to the 
symmetry of the Gaussian field ho- We will ignore the 
last term since our analysis is restricted to first order 
in A/4za All that remains is the two-point correlation, 
which with the help of eq. ([2]) is seen to be 

(vo(h)vo(f 2 )) = (ho(fi)h (f 2 )) 

= ^±A(fc) 2 cos(£- (h -f a )). (30) 



We will now plug our intermediate results, eqs. (|28[) and 
(|30)) . back into eq. (|27p . For convenience, we will set r = 
0, which we are allowed to do thanks to the homogeneity 
of v. We find 

a = J2 IMK) 2 1 1 d 2 f x d 2 f 2 n- 2 (A V t)-\h ■ f 2 ) 



cos(/c • (fi - f 2 )). 



(31) 



Note that based on eq. (|28j) we should have put {x\ — 
iyi){x 2 + iy 2 ) instead of (f\ ■ f 2 ); the latter is merely 
the real part of the former. However, since we already 
know that the final answer is real (since a = (|u z | 2 )), we 
can conclude that the imaginary part would not give a 
contribution. 

After performing the integrals in eq. ([3~Tj) we get the re- 
sult given below. The three-point correlations /3, 7 and S 
give rise to six-dimensional integrals involving four-point 
correlations (which are first order in X/4i/). These corre- 
lations can be factorized into two two-point correlations 
by Wick's theorem, resulting in a sum over two wave vec- 
tors k\ and k 2 , as opposed to the one we had in the case 
of a. 

All the relevant correlations are 



4 -2k z vt 



(32a) 
(32b) 



4z/ W 4 



xe 



-2(k 1 1 +kl+k 1 -k 2 )vt 



(32c) 



jr, E E \m?m?^ti&&i + k 2 f 



4f W 4 

Aii k 



xe 



A EEi^) 2 ^ 2 ) 21 



4^ w 4 

Aii &2 



32 

2{k'l+k'^+k 1 -k2)vt (32d) 

k\k\(k\ + k 2 ) 



32 



+ ((fci+fc 2 ) 4 -fc 1 -fc 2 1 )(fci-fc 2 ) 

2(fcf+fc|+fcl-fel)i/* ($2c) 



xe 



6 



< 




duces nontrivial functions. In this case, eq. (|17|) reads 



0.05 



t (x vkj 

(a) 




t (x VkJ) 

(b) 



Figure 2: The imbalance between maxima and minima An 
of h(f,t), where h obeys the KPZ equation (with A/4f — 
0.1), as a function of time. At t = 0, /i(r) was taken to 
be a Gaussian field with (a) a Gaussian spectrum A(k) oc 
exp(— fc 2 /(4fco)); (b) a ring spectrum A(k) oc 5(k — ko). Shown 
are our theoretical perturbative result (eq. {T7}) and data 
from simulations. 



For a continuous spectrum, the sums can be replaced by 
integrals. 

We see that the parameters depend on the spectrum of 
ho in a nontrivial way. Especially the presence of k\ ■ k% 
(which is also present in terms like (k\ +/C2) 2 ) in the rela- 
tions for 0, 7 and S complicates matters, as it introduces 
a dependence on the angle between k\ and k%. An exact 
analytical evaluation is therefore only realizable for a few 
spectra of a convenient form. Even for the so-called ring 
spectrum, with A(k) oc 6(k — ko), arguably the simplest 
spectrum one can have, the angular dependence intro- 



An= F 7\ (2 + r)/ (2r) - 5rh (2r) 

4y 9 V ir t V (33) 

+ (2 + r + 6T 2 ) Fi(2;r 2 ) N ), 



where r = k^vt] I and l\ arc modified Bessel functions 
of the first kind and qFi is the confluent hypergeomctric 
function. Recall that we set Ko = (/i 2 ,) = 1 for conve- 
nience; for the general case, a factor of y/K~o needs to be 
added. 

Another, more elegant case in which an exact evalu- 
ation of eq. p7|) is possible is the Gaussian spectrum 
A(k) oc exp(-fc 2 /(4fco)), for which 



An 



A 



64t 3 (1 + 4t) 7 /2 



4z/ V3^(1 + 2t) 3 (1 + 6t) 4 ' 



(34) 



where again r = k^vt and a factor of y/Ko needs to be 
added for our result to apply in general. 

Going back to the general case of an unspecified power 
spectrum, it is convenient to expand An in t. The result 
is 



An 



Al 
4^9 



6 2K 2 K 6 



"\K 2 



0(t 3 



(35) 



for all Kq- One may note that for a Gaussian spectrum, 
there is no quadratic order in eq. (|34p . which is confirmed 
by the above formula, since 1K<xK§— 3A"| = in this case. 

The analytical results for An above are compared to 
results from numerical simulations (with Ko = 1 and 
A/4z/ = 0.1) in figs. 2(a) and 2(b) The general method 



is the same as outlined in ref. |7j|. We start with a Gaus- 
sian field ho defined on a finite square grid with periodic 
boundary conditions. We then transform to vo and use 
the alternating direction implicit (ADI) method to simu- 
late the heat equation, collecting statistics on the maxima 
and minima at every time step. The results are averaged 
over for tens of thousands of /lo's, each with the same 
spectrum but random phases. 

In general, if a field evolves under a nonlinear equa- 
tion for a long time, the non-Gaussianity can become 
large, even when the perturbation is small, because it 
will add up over time. Thus we may expect a breakdown 
of our predictions after some time, as in fig. 2(b) How- 
ever, the KPZ equation has a special mapping to a dif- 
fusion equation (eq. (|22p). and this implies that the non- 
Gaussian perturbations never build up. Eq. (|26|) shows 
that the nonlinear correction diffuses outward but does 
not grow over time. Therefore, for the KPZ equation, 
our approximations should remain accurate for arbitrar- 
ily long times. This is indeed what we see in fig. 2(a) 
where h(t = 0) is Gaussian field with a Gaussian spectral 
function. 



In fig. 2(b)| however there is a breakdown for the ring 
spectrum. This spectrum is special because it has zero 
weight at k = 0. This implies that the leading Gaussian 



7 



term in eq. (|26[) is suppressed exponentially, decaying as 
cxp(— k^vt) (see eq. (|2"Tj)). Thus after a long time, the 
second term dominates, and our approximation that v 
is close to a Gaussian no longer holds. Whenever the 
spectral function has a weight at k = (as in fig. 2(a) I, 
the approximation works for a longer time. 



III. CONCLUSIONS 

We have found a general perturbative formula, eq. (1171) , 
for determining the imbalance between maxima and min- 
ima of an isotropic random field that is almost Gaussian. 
It allows one to attack the reverse problem, namely, to 
determine the size of the phenomenon that causes the 
non-Gaussianity, by measuring the relative densities of 
maxima and minima. In the case of the deterministic 
KPZ equation for instance, the imbalance can reveal the 
size of the nonlinear parameter A relative to the diffusion 
coefficient v. 

In ref. Q , we investigated the imbalance between max- 
ima and minima as a result of non-Gaussianity. Although 
we arrived at an exact result, it applied only to the spe- 
cial case of a local perturbation, i.e. for a field given by 
h(f) = FMh{H{f)) where H is a Gaussian field and Fnl 
any (nonlinear) function. The result in the present study, 
although perturbative, also accommodates nonlocal per- 
turbations, provided that the resulting field is still homo- 
geneous and isotropic. 

For local perturbations, we found that the size of the 
imbalance is exponentially small in the size of the per- 
turbation Q. Nonlocal perturbations however allow for 
a power-law relation. This is apparent in eq. Q35p , which 
shows that the KPZ equation can cause an imbalance 
that grows quadratically with time. As a result, the den- 
sities of maxima and minima can prove to be a sensitive 
test to not only detect non-Gaussianity, but also to dis- 
tinguish local from nonlocal perturbations that induce 
non-Gaussian statistics. 
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Appendix A: Higher order cumulants 



Consider the equation 

h n = 'S ' A nm + £ S ' B n pqhphq, (Al) 



with the initial condition 



h n (0) = H n , 



(A2) 



where the H n 's are a set of variables with a joint Gaus- 
sian distribution. These coupled differential equations 
are a simple model of a nonlinearity, with the lowest or- 
der (quadratic), and they also include the KPZ equation 
as a special case, if it is discretized. This differential 
equation illustrates the general principle that cumulants 
of a high order are very small if the nonlinear term in 
the differential equation is small - unless one waits long 
enough for these cumulants to build up. 

For this family of equations the precise result is that, 
after a finite period of time, the fc-th order cumulants of 
any of the h n 's are of order at most e k ~ 2 if k > 2 (for 
k = 1 or k = 2 they are bounded). 

There are two steps in the proof: first, we find how h n 
depends on the initial conditions, and show that it has 
the form of a power series in e. The result is that 



h n (t)=Fi°H{H j })+ S FU({H j }) 
+ s 2 F^({H J }) + ... 



where Fn°^ is a linear function, F^ 1 



(A3) 



is quadratic, etc. 
So the dependence of a given term on the Hj 's is poly- 
nomial; the dependence on t is all in the coefficients of 
these polynomials. 

In other words, h n can be expressed in the form of a 
nonlinear function of a Gaussian, the same type of func- 
tion whose cumulants we calculated in We will see 
that many of the cumulants vanish - this is the second 
step of the proof. We calculate the cumulants, 



C k (h ni ,...,h nk )=J2e r E C ^ F n? } 

r=0 r 1: r 2 ,...,r k 



(A4) 

All the terms up to order r = k — 3 vanish, so that the 
remaining terms are of order e k ~ 2 or smaller. This is 
a consequence of a general theorem: a cumulant of k 
polynomials in Gaussian variables is zero if 



k > 1 



d 
2' 



(A5) 



where d is the sum of the degrees of the polynomials. 
In Ck{Fni\ ■ ■ ■ , F„ h ) the sum of the degrees is d = 
J2i n + 1 = r + k. If r < k - 3, then Eq. (|A5]) follows, 
so the cumulant vanishes. 



In this section it is demonstrated that, for an initially 
Gaussian field evolving according to a diffusion equation 
with a perturbative nonlinear term, the cumulants be- 
come smaller as the order increases (i.e. they are of higher 
order in the perturbation). 



1. Power series solution 

Expand h n (t) = e r hn\t) and substitute it into 
eq. (|A1[) . and then match the coefficients of e r . This 
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gives the relation 

m 

r-1 (A6) 

=EEM nl ( t )^ 1 (*)' 

p,q T\= 

Here, everything depending on h^ r > is on the left-hand 
side; everything on the right-hand side depends on earlier 
terms in the series, with 7'i < r. This means that 
one can solve the equations recursively: first hnd the 
h's up to r! = r — 1, then substitute it into the right- 
hand side of the equation and then solve for h^ r \ which 
is straightforward because it is a linear equation with 
a source. We only need to know the initial conditions, 
which are 

= H n ; h n r] = for r > 1. (A7) 
The solutions to the equations are given as follows: 

h^(t)=^[e M ] nm H m , (A8) 



where e At is the exponential of the matrix At, which is 
just a set of functions of t. 

These functions are all polynomials in the -ff/s. First, 

h n is obviously linear. Entering r = 1 in eq. (|A9[) shows 
that hn is the sum and integral of hp^h q \ which is 
thus quadratic in the ifj's. Now we can find the general 
dependence inductively: assume that it has already been 
shown that hn is a degree r\ + 1 polynomial in the Hj's 

for n < r. Then ^V)^" 1 "^^) is of degree r + 1, 
and thus h^ r ' is as well. 



2. Vanishing cumulants 

We will calculate the cumulants of polynomials in the 
iij's by reducing them to cumulants of the -ffj's them- 
selves, which are Gaussian. A helpful identity for this 
expresses C{xy, z\, . . . , z q ) where x, y, z% are any random 
variables in terms of simpler cumulants. The identity is 

C(xy, Z!,...,z q )= C(x, y,z 1) ...,Zq) 

+ C(x,z s )C(y,z T ). ( Al °) 

SUT={l,...,q} 

The sum is over all ways of partitioning the indices of the 
z's into two sets S and T. The symbol z$ is short for the 
list of all the z's corresponding to the indices S. 



Here is an example: 

C{xy, u, v) = C(x, y, u, v) + C{x)C{y, u, v) 

+ C(x,u)C(y,v) + C(x,v)C(y,u) (All) 
+ C(x,u,v)C{y). 

A proof of this relation can be obtained using induc- 
tion. First note that it is trivially true for q = 0, since 
C(x,y) = (xy) — (x)(y). Now we assume the relation to 
hold for all q' < q. Consider the identity (see e.g. Q or 
0) 

(x 1 ...x n } = ^C(x Sl )C(x S2 )...C(xsJ, (A12) 

where the sum is taken over all the ways in which the 
set {1, . . . , n} can be partitioned into disjoint subsets Si- 
If we apply this identity to the set {x, y,Z\,..., z q } and 
group together the terms for which x and y are in the 
same subset or in different ones, we find 

(xyzi ...z n ) 

= ^2 c { x ^y^ z u)c{z Vl ) . . .c{z Vm ) 

+ Y, C(x,z s )C(y,z T )C(z Vl )...C(zv m ) 

S,T,{V % } 

= Y C(z Vl )...C(z v j[c(x,y,zu) (A13) 

U,{Vi] 

+ C(x,z s )C(y,z T ) 

SUT=U 

We can also choose to expand (xyzi . . . z n ) while treating 
xy as a single variable, which results in 

(xyzi . . . z n ) = Y c ( x y> z u)C{z Vl ) ...C(z Vm ) 
u,{Vi} 

(A14) 

The two decompositions into cumulants should be equal. 
By induction, we can pose 

C(xy,zu) = C(x,y,zu)+ ^ C(x, z s )C(y, z T ) (A15) 

SUT=U 

for all U ^ {1, . . . , g}. It then easily follows that the 
relation must also hold for U = {1, . . . , q}. 

We will use this identity to prove that if p\ , . . . pk are 
degree d\,. . .dk polynomials in Gaussian variables and 
d = J2idi, then Ck(pi, ■ ■ ■ ,Pk) vanishes if cq. (|A5|) is 
satisfied. We shall first demonstrate the procedure us- 
ing a simple example: C(H 2 ,H 2 , H, H, H) where H is a 
Gaussian variable. We will reduce this to cumulants of 
H by using eq. (|A10[) ; that will mean we have to apply 
the identity twice to split up both H 2, s. After the first 
time, we have a sum featuring one term with a single cu- 
mulant, C(H 2 ,H, H, H, H, H), while the other terms are 
products of two cumulants. Furthermore, there is only 
one H 2 left in each term. After applying eq. (|A10|) a 
second time, we are left with products of at most three 
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cumulants. Since there are 7 H's distributed among these 
cumulants, at least one of the cumulants in each product 
is of at least third order, and hence zero because the H's 
are Gaussian. Hence C(H 2 , H 2 , H, H, H) = 0. 

In general, we first use the multilinear property 
of the cumulant function (i.e. C(x + y, z, w, . . .) = 
C(x, 2, w, . . .) + C(y, z, w, . . .)) to reduce each of the vari- 
ables to one term (which is a product of some of the H's). 
It takes d—k applications of eq. (|A10I) to split all the vari- 
ables up into individual H's, because it takes di — 1 steps 
to factor the z-th variable, for a total of JT di — l = d — k 



steps. Since each application of eq. (|A10[) adds at most 
one cumulant to each term, in the end each term has 
at most d — k + 1 factors of C. This is less than | by 
eq. (|A5[) . But there are a total of d variables H's that are 
split among them. Hence one of the factors is a third- 
order cumulant or higher, which means that it has to be 
zero. 

Now this result can be combined with eq. (|A3|) to prove 
that the fc-th order cumulants of the h n 's are of order 
e k ~ 2 , as we showed above. 
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